Planform selection in two-layer Benard-Marangoni convection 

A. Engel* a and J. B. Swift 6 

a Institut fur Theoretische Physik, Otto-von-Guericke Universitat, Postfach 4120, D-39016 Magdeburg, Germany 
Center for Nonlinear Dynamics, University of Texas at Austin, Austin, TX 78712, USA 

Benard-Marangoni convection in a system of two superimposed liquids is investigated theo- 
retically. Extending previous studies the complete hydrodynamics of both layers is treated and 
buoyancy is consistently taken into account. The planform selection problem between rolls, squares 
and hexagons is investigated by explicitly calculating the coefficients of an appropriate amplitude 
equation from the parameters of the fluids. The results are compared with recent experiments on 
two-layer systems in which squares at onset have been reported. 
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The hexagonal convection cells discovered by Benard in his famous experiments on thin oil layers heated from below 
ID have become the trademark of pattern formation in hydrodynamic systems driven slightly out of equilibrium (see 
e.g. P|). The hundred years of research devoted to this system have revealed several important insights but also wit- 
nessed several misconceptions. Rayleigh's original theoretical description |3| focusing on buoyancy-driven convection, 
£S| \ though indicating a possible instability mechanism, failed to produce a threshold compatible with experiment. Not 
^ ■ until forty years later was it realized that the temperature dependence of the surface tension is the crucial driving 
force in thin layers B. The corresponding linear stability analysis jq] gave stability thresholds consistent with the 
experimental findings, moreover, a subsequent weakly non-linear analysis |6j,[7| produced theoretical support for a 
sub-critical transition to a hexagonal flow pattern M. 

Quite naturally the first theoretical investigations were performed using simplified models of the experimental 
Q\ • situation. The initial assumption of a flat surface of the liquid was soon relaxed by Scriven and Sternling and 
^\ | Smith |lffl who were able to show that surface deflections give rise to an additional instability appearing at very long 
' | ■ wavelengths. It was only very recently that this instability was unambiguously demonstrated in an experiment pd| ] 
where it manifests itself as a distortion of the layer thickness with a characteristic length which is of the order of the 
lateral extension of the fluid layer. Being observable only in very shallow liquid layers, the instability usually results 
~j—{ ' in the formation of dry spots. 

Another common simplification is the restriction of the instability mechanism to either buoyancy or thermocapil- 
Q-r larity [0,0,[l8[ , although there seem to be rather few experiments M|||]T[] which have been performed in parameter 
"> regions with the ratio between the Rayleigh and the Marangoni number being sufficiently different from unity. Also, 

most investigations focussed on a single layer model in which a lower liquid layer is in contact with a gaseous upper 
layer and only the hydrodynamics of the liquid are treated. The convection in the gas layer is neglected and the heat 
exchange between the layers is often modeled in a phenomenological way using a Biot number, see e.g. p5f. Even 
if a genuine two layer model is considered the viscous stresses and the pressure variations in the gaseous layer are 
neglected in order to keep the analysis simple fll3[ . 

On the other hand it has been known for some time fljjlll that a system of two superimposed liquids displays 
a much richer behaviour than the single layer models. In particular the Marangoni instability can be induced by 
heating from above such that buoyancy and thermocapillarity compete rather than enhance each other, a situation 
which in single layer systems can only be realized using the rare case of liquids with anomalous thermocapillary effect 
in which the surface tension increases with increasing temperature U%- Many additional features such as oscillatory 
instabilities |18| or transitions from up- to down-hexagons may be found in systems with two liquid layers. The rich 
variety of phenomena occurring in the theoretical analysis of the two-layer liquid systems results in part from their 
huge parameter space. A single layer system is characterized by just three dimensionless parameters; namely, the 
Rayleigh number, the Marangoni number and the Prandtl number. The latter is irrelevant in the linear analysis and 
the first two are both proportional to the temperature difference across the layer. Two layer systems on the other 
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hand may easily need ten or more dimensionless parameters for a complete specification. These numbers include the 
ratios of the hydrodynamic parameters of the participating liquids. 

For a long time Marangoni convection in two-liquid-layer systems was an interesting theoretical problem but too 
difficult to handle experimentally. Already Zeren and Reynolds |12f| tried to experimentally realize the instability by 
heating from above which came out of their theoretical analysis but failed. Very recently, however, experiments where 
performed in which the Marangoni instability in 1-2 mm thick superimposed layers of immiscible liquids was observed 
pl| , p2| . In particular an instability by heating from above and square patterns at the onset were reported. 

In the present paper we will investigate theoretically Benard-Marangoni convection in a system of two liquid layers. 
Building on the linear stability theory developed in [G3[ we perform a weakly non-linear analysis in order to solve 
the planform selection problem slightly above the linear stability threshold. To this end the competition between 
rolls, squares and hexagons will be discussed. Only perfect patterns will be considered leaving the question of 
weakly modulated patterns for future investigations. We will consistently include buoyancy effects and treat the full 
hydrodynamics of both liquids, generalizing in this way various previous treatments [p 13 ; 13,E4 E6[. However, we will 



assume a flat interface between the two liquids. As will become clear below, interface distortions are crucial for the 
long wavelength instability resulting in dry spots but can be safely neglected when dealing with the finite wavelength 
instability resulting in cellular patterns. 

The paper is organized as follows. In section 2 the basic equations are collected and transformed into a form 
suitable for the weakly non- linear analysis. Then the perturbation scheme is set up and the necessary computational 
steps are listed. Section 3 deals with the first order of the perturbation theory which is nothing but the linear 
stability analysis. In section 4 the main steps of the nonlinear analysis are outlined. The solution of the second order 
problem is relegated to appendix C and the solvability condition in third order is then formulated to derive the desired 
amplitude equation characterizing the planform selection problem. Section 5 discusses the results obtained for some 
experimentally relevant combinations of liquids. Finally section 6 contains a discussion of the results together with a 
comparison with experimental findings. 

II. BASIC EQUATIONS 

We investigate a system of two layers of immiscible and incompressible liquids of thickness ffi' with densities p^\ 
kinematic viscosities v^\ coefficients of volume expansion cS l \ heat diffusivities x > an d thermal conductivities k> 1 > 
where the superscript i = 1 (2) denotes the lower (upper) fluid (see fig-lO). The system is bounded in the vertical 
direction by two solid, perfectly heat conducting walls with fixed temperatures T b and T* and is infinite in the 
horizontal directions. The interface between the two fluids is assumed to be flat and to lie in the x-y-plane of the 
coordinate system. 
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FIG. 1. Sketch of the system under consideration. One liquid layer is superposed on another between two horizontally 
infinite, perfectly heat conducting plates. The interface between the liquids is assumed to be flat. Convection arises due to 
buoyancy and the temperature dependence of the surface tension. 

The hydrodynamics of the two liquids will be described within the Boussinesq approximation, i.e. we assume that 
all parameters are independent of the temperature, except for the densities p^' and the interface tension a. More 
precisely we use p {i) (T) = p® (T b )(l-a^ (T-T b )) and V ± er = da/dT V ±T with constant a (i) and da/dT. Neglecting 
heat production due to viscosity, the basic equations describing the system are the continuity equations 

Vt> w =0 , (2.1) 



the Navier-Stokes equations 

dtv® + (v^V)v^ = --J^Vp w - ,g(l - a w (T« - T b ))e z + z/«Ai>« , (2.2) 

and the equations of heat conduction 



p(i) 



t r« + (« (i) v)r« = x «atW . (2.3) 

Here e z denotes the unit vector in the vertical direction and g is the acceleration due to gravity. 
These equations are completed by the boundary conditions 

„(i) = o , T« =T b at z = -/i (1) , (2.4) 

and 

„(2) = f T m = T t at JS = ft (2) ; (25) 

at bottom and top respectively and 

w (i) =t; (2) j T (i) =T (2) ; K V-)d z T<U = K^d z T^ 

(a^-a^)e z ]^ = -^ ± T , v™ = v™ = , at z = , (2.6) 

expressing the continuity of the velocities, temperatures and heat fluxes respectively as well as the balance of tangential 
stresses at the interface. The a^ 1 denote the stress tensors in the liquids and the subscript _L describes the projection 
to the x-y-plane. In accordance with our assumption of a flat interface between the liquids the condition for the 
continuity of the normal stress at the interface is replaced by the requirement that the perpendicular components of 



the flow velocities must vanish. This is expressed by the last equation in (2.6) . 

Introducing ft.' 1 ), (ft^) 2 /x > X A an d p^ v^- 1 ' X / (h ) 2 as units for length, time, velocity, and pressure 
respectively we find for the velocities v — (u, v, w) (V = (U, V, W)) and the appropriately normalized deviations (0) 
of the temperatures from their static profiles in the lower (upper) liquid the equations: 

— (d t v + (v\7)v) = -X7p+R6e z +Av (2.7) 

Pr 

d t 6+{vV)6 = w + A9 (2.8) 

— (d t V + (W)V) = -VP + aR9e z + vAV (2.9) 

Pr 

d t @ + (V V)0 = - W + xA9 , (2.10) 

K 

where the pressure fields p and P in the lower and the upper liquid differ from p^ 1 ' and p( 2 > respectively only by trivial 
contributions stemming from the buoyancy terms. The boundary conditions acquire the form 

v = Q , 6> = at z = -l , (2.11) 

V = , 6 = at z = a , (2.12) 

and 

v_L = V_L,w = W = 0,9 = Q,d z 6 = Kd z @,d z w-r]dzW = MA±e, &tz = 0, (2.13) 

where in the last equation the continuity equation was used. Moreover the following parameters have been introduced: 

h {2) a (2) „(2) p {2) K (2) (2) 

a= ^' a= o^)' l/= ^)' r ^^' K= ^)' X= ^)' (2 - W) 

as well as the Prandtl-number Pr = u^> /x ■> the Rayleigh-number 

R= V W X W a + K {T " T) ' (2 - 15) 



and the Marangoni-number 

M = -% .Jfl (1) —(T»-Tt) . (2.16) 

For the Rayleigh- and Marangoni-number we have chosen the standard expressions corresponding to the lower liquid. 
The respective numbers for the upper liquid are then given by 

R m = ^ R and M (2 ) = ^!_ M (217) 

vxk XV K 

respectively. 

The ratio between the Rayleigh and Marangoni numbers determines whether the occurring instability is predom- 
inantly driven by buoyancy or by surface tension. Experimentally both parameters are varied simultaneously since 
they are both proportional to the temperature difference T b — T* . We will therefore replace R by cM with the 
temperature independent constant 

specifying the experimental setup. In this way both buoyancy and surface tension are included in a consistent way. 
We assume that da/dT < as is the case for most systems of two liquids such that c > 0. Note that both the situation 
of heating from below and heating from above are described with the latter case corresponding to M < 0. 

The set of equations may be simplified by standard manipulations. Taking twice the curl of the Navier-Stokes 
equations, using the continuity equations, and projecting onto e z we get the following basic set of equations for the 
^-components of the velocities and the temperature fields: 

A 2 w + cMA ± 9 = -j- [d t Aw - d z (V ± {v\7)v ± ) + A ± {vV)w] (2.19) 

Pr 

w + AQ = dtO + («V)0 (2.20) 

vA 2 W + acMA^Q = — \d t AW - <9 Z (V AW)V ±) + A ± (VV)W] (2.21) 

Pr 

-w + xAe = d t e + (vv)q (2.22) 

K 

together with the boundary conditions 

(2.23) 
v d 2 z W = MA ± 6 at z = (2.24) 

(2.25) 

In order to investigate the planform selection problem we will derive third order amplitude equations for the slow time 
variation of the amplitudes of different unstable modes. Similar to the case of the Rayleigh-Benard instability B] the 
no-slip boundary conditions at top and bottom suppress the vertical vorticity, i.e. (V x v)e z = (V x V)e z = 0, and 
therefore w e do n o t exp ect problems due to a coupling to a slowly varying mean flow |28| up to this order. From the 
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solution of ( 2.19 )-( 2.22 ) we hence obtain w, 9, W and 8. Using the continuity equations and the absence of vertical 
vorticity allows to determine u, V and U, V and finally the pressure fields follow from the Navier-Stokes equations. 
It is convenient to write the above equations in the form 

L<p = T(<p)+N(ip,ip) (2.26) 

with the state vector 

<P=\w\ ' ( 2 - 27 ) 



and the linear operator L defined by 



L = 
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(2.28) 
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and the boundary condit ions ( 2.2S )-( 2.25 ). T(tp) d enote s the time dependent terms and ftf(<p,(p) describes the 
quadratic nonlinearity in (2.1£)-(2.22). We will solve ( 2.26D perturbatively using the ansatze 

(2.29) 

(2.30) 

(2.31) 

with a small parameter e. In the case of a static instability we have to — whereas for an oscillatory instability w^O 
gives the frequency of oscillation of the unstable mode. Using the perturbation expansion specified above we consider 
a situation slightly above the threshold M c o f the li near i nstab ility, where the amplitude of th e unstable modes can 
still be considered to be small. Plugging ( 2.29 )-( 2.31 ) into ( 2.26 ), taking into account that ( 2.30 ) implies an expansion 

L = L + sL 1 +e 2 L 2 + ... (2.32) 

for the linear operator and matching powers of e the non-linear problem transforms into a sequence of linear equations 
of the form 



L <Po = 

Lopi = -Lupo + Af{ip , ipo) 

L <P2 = -L 2 <Po - Lupi +T(<p ) +J\f(ipi,ip ) +Af(ip ,fi) 



(2.33) 

(2.34) 
(2.35) 



The first line is just the linear stability problem. The condition for non-trivial solutions ipo of this equation makes 
Lq singular and yields the critical value M c of the bifurcation parameter M. From the translation invariance in the 
z-y-plane we know that ipo is of the form 



^o = Vo(-s) exp(ifcr — iuoi) 



(2.36) 



where r — (x,y) and k — (k x ,k y ) are two-dimensional vectors. There is a critical value M c (k) of the bifurcation 
parameter for all values of \k\ = k and minimizing M c (k) in k gives the wavenumber k c of the first unstable mode 
together with the critical Marangoni number M c = M c (k c ). 



The remaining equations in the hierarchy starting with (2.34) all involve the very same singular operator Lq but are 
mhomogeneous. Consequently the perturbation expansion makes sense only if the inhomogeneities are perpendicular 
to the zero eigenfunction of the adjoint operator Lt of Lq. 
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FIG. 2. Relative orientation of the two-dimensional wave vectors appearing in the ansatz (2.37). The two triads fei, &2, &3 
and ki, fcs, &6 of wave vectors wi th fcs perpendicular to fei allow to describe rolls as well as squares and hexagons by different 
values for the amplitudes A n in (2.37). 



In order to address the planform selection problem within the perturbation approach sketched above the form of 
<po must be sufficiently general and in particular must include the different planforms observed in the experiment. We 
will discuss the planform selection problem only for the case of the static instability leaving the investigation of the 
oscillatory instability to future work. It is then sufficient to use for ipo the form 



ipo = ipo(z) 



J2Mr)e lk " r 



(2.37) 



with the six two-dimensional vectors k n obeying \k n \ — k c and fci + k 2 + &3 = 0, k 4 + k$ + k$ = 0, as well as k\k§ = 
(see fig.0). Depending on the values of the amplitudes A n this form describes rolls (e.g. A\ — A, A n = for all n > 1), 
squares (e.g. A\ — A 5 = A, A n = else) and hexagons (e.g. A\ = A 2 = As = A, A n = for n > 3). 

Using this form we find from the solvability conditions of (2.34) and (2.35) an equation describing the time evolution 
of the scaled amplitudes A n — sA n . As is well known |2j the general form of this amplitude equation already follows 
from the symmetries of the problem. For the present situation it is given by 



d t A 1 = eAi + 7A2A3 - \Ai 
with the super-criticality parameter 



9h(\A 2 f + \A 3 \ 2 ) + g t (\A 4 \ 2 + |i 6 | 2 ) + g n \A 5 \ 2 A, 



M - M r 



M c 



(2.38) 



(2.39) 



Similar equations for the other amplitudes follow from permutation and complex conjugation. The terms included in 
these equations are the only ones up to third order which are invariant under the transformation A n 1— ► A n exp(ifc„T*o) 
corresponding to a translation by r in the x-y-plane. Moreover due to the isotropy in the x-y-plane the coupling 
coefficients between the different terms in (2.37) may only depend on the angle between the corresponding wave 
vectors. 



A well known linear stability analysis of the various fix points of ( 2.38 ) yields the stability regions of the different 
planforms as functions of the parameters e,j,gh,gt,gn p7| |. The remaining problem is hence to use the perturbation 
expansion described above to express these coefficients of the amplitude equation in terms of the hydrodynamic 
parameters of the problem. To this end the following well-known program has to be carried through: 



• Calculate M c (k) from the linear problem and determine k c = argminM c (fc) and M c 

• Determine the adjoint operator Lq of Lq and its zero eigenfunction <pq. 



M c (k c 



Calculate the inhomogeneity of the 0(e 2 )-equation ( 2.34 ) and apply the solvability condition to this order. 

Solve the 0(e 2 )-equation ( p. 34 ) to determine <p\. 

Calculate the inhomogeneity of the 0(e 3 )-equation ( |2.35| ) (only terms proportional to expi(fcir) are necessary) 



• Combine the solvability conditions at order 0(e 2 ) and 0(e 3 ) to derive ( 2.38 ) and extract the expressions for the 
parameters 7, g h ,g t ,g n . 



III. THE LINEAR PROBLEM 



We first solve the 0(e) problem ( 2.33| ), which is equivalent to the linear stability analysis. Putting 
¥>o = fo{ z ) exp(ifcr — imt) and using the ansatze 



we find 



w (z),6o(z) ~ exp(Az) 



Wo(z),e (z)~exp(A*) 



(3.1) 



(A 2 - fc c 2 )(A 2 - k\ + ^)(A 2 - k 2 + »w) = -cMk 2 (A 2 - k 2 c )(A 2 - k 2 + ^-)(A 2 - k 2 -) - -—cMk 2 c 
Pr vPr x V ^X 



(3.2) 



We therefore obtain six different values for Ai and Aj. It is convenient to define Xi — A^q for i 
write 



12 and to 



wo(z) = ^2w 0i e * 2 
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(3.3) 

(3.4) 



The boundary conditions (2.23)-(2.25) give then rise to a homogeneous system of linear equations for the 12 unknowns 
Woi- In order to get a non-trivial solution the determinant of the coefficient matrix A must vanish. The conditions 
for the real and the imaginary part of det„4 yield the desired functions M c (fc;par) and w c (fc,par) where par = 
(a, a, k, x, v, rj, c, Pr) stands for the vector of parameters in the problem. 
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FIG. 3. Dispersion relation M c (k) as resulting from the linear stability analysis for the hydrodynamic parameters of setup 
2 listed in appendix A. The system shows an instability when heated from below as well as one when heated from above. 

A typical result for a static instability is shown in fig.0 displaying the dispersion curve resulting from the numerical 
analysis of det.4 = for u = using the parameters of setup 2 listed in appendix A. As can be seen from the figure in 
this system one may have an instability by heating from below (M > 0) as well as when heating from above (M < 0) . 
In fig|| the results of the present approach for the setups 1 and 5 of appendix A are compared with those resulting 



from the full linear stability analysis including surface deflections as considered in 23 1 . As is clearly seen in the region 
of the pattern forming instability k = 1...3.5 the two curves are almost identical with differences showing up only for 
small wave numbers k ^C 1. Within the linear theory the surface deflections for unstable modes corresponding to the 
planform selection problem may therefore safely be neglected and we expect that this is also a good approximation 
for the weakly non-linear regime. 
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FIG. 4. Comparison of the dispersion relation M c (k) as resulting from the present analysis assuming a flat interface between 
the liquids (thin full lines) with the results of the complete linear analysis of f23| including surface deflections (dashed lines). 
The lower right curves are for setup 1, the others for setup 5 specified in appendix A. 



Having obtained the dispersion relation we calculate k c by minimizing M c {k) and determine the critical Marangoni 
and Rayleigh numbers of both fluids as well as the temperature difference across both layers at the instability. The 
results for the setups under consideration are summarized in the upper part of table 1. 

From all the parameters of the system the depth ratio a is the only one which may be easily varied in the experiments. 
For the parameters of setup 3 and a total depth of 4.5 mm we have calculated the critical Marangoni number and the 
critical wave number modulus as a function of the thickness hS l > of the bottom layer restricting ourselves to the case 
of heating from below but including the possibility of an oscillatory instability. The results are displayed in figjq. For 
values of ft/ 1 -* between 1.5 and 2.5 an oscillatory instability precedes the static one which would occur at unusually 
large Marangoni numbers only. A similar oscillatory instability was also found for a two-layer system in which the 
Marangoni effect was neglected and pure buoyancy-driven convection was considered, and an intuitive interpretation 
as a periodic change between viscous and thermal coupling of the flow fields at the interface was given |18|. The 
oscillatory instability was also detected in the experiment using setup 3 with M 1 ' = 1.8mm and the experimental 
values for the critical Marangoni number and the wavelength of the oscillatory mode are in good agreement with the 
theory 62j. 
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FIG. 5. Critical Marangoni number for heating from below a system with parameters as specified in setup 3 of appendix 
A and total dep th 4.5 mm as a function of the bottom layer thickness hr*. Note that both M and k are scaled with hr- 1 ' (cf. 
(Eul and (Eja)). 



Knowing the critical value of M we can now also determine the coefficients of the eigenvector corresponding to the 
zero eigenvalue. This fixes the functions Wq(z), 0q (z), Wq(z) and &o(z) up to an overall constant and completes the 
determination of ipo. 

Finally we have to consider the adjoint problem and to calculate its zero eigenfunction tpo where we again restrict 
ourselves to the stationary instability. The adjoint operator L + is determined in appendix B. The calculation of its 
eigenfunction to the eigenvalue zero is very similar to the determination of ipo described above. We find that it is of 
the form ipo exp(ik n r) where the components of ipo may be written as 



w (z) = ^2w Ql e iJ 

i=l 
12 

W (z)=J2™°i eXiZ 



w Ch c A lZ 



0o(z) = cMk^ x2 _ k2 

i—l l c 



12 



acMk 2 c ^ 



w 0i A, 



X t^ X *- k c 



(3.5) 
(3.6) 



with the same parameters A.; as determined by ( p.2|) with u> — 0. The boundary conditions give again rise to a 12 x 12 
system of linear homogeneous equations for the coefficients woi- As before the condition for a non-trivial solution 
is a vanishing determinant of the corresponding matrix. Note, however, that there is now no parameter to adjust! 
The deviation of the smallest eigenvalue of the matrix found in the numerical calculation from zero gives therefore a 
valuable hint on the accuracy of the numerical procedure employed. 



IV. THE NONLINEAR ANALYSIS 



The solution of the planform selection problem requires the treatment of the nonlinear interaction between different 
unstable modes. To include nonlinear terms up to the third order in the amplitudes A n introduced in (2.37) we have 



first to solve (2.34). The general procedure is standard, some intermediate steps are sketched in appendix C. Using 
this solution we are in the position to calculate the terms appearing on the right hand side of ( p. 35 ). We do not 
have to solve this equation, but only need to know the solvability condition at this order. Due to the x-y-integrals in 
(B9) and the r-dependence of ipo only terms proportional to cxp(±ife„r) give rise to non-trivial contributions to the 
solvability condition. In fact it is su fficien t to focus on terms proportional to exp(ifcir) since these finally give r ise to 
an amplitude equation of the form ( 2.38J ) for A\. Equivalent equations for the other amplitudes of the ansatz ( 2.37 ) 
follow then from permutation and complex conjugation. 



In order to collect the relevant terms we first realize that there are contributions 
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(4.1) 



originating from the terms —L 2 ip , —L\<p\, and T(tp ) respectively in (2.35). Here Q\ and Oi denote the solutions 
obtained in the last section for the resonant term. 



The contributions proportional to exp(ifcir) from the last two terms in ( 2.35 ) arise from combinations between 
ip ~ exp(iqr) and <px ~ exp(ipr) with q + p = k±. From the continuity equation, \7v — 0, and the absence of vertical 
vorticity, (V x v)e z = 0, we find 



iar ^Q q 

VOX = e q -J OzWo 
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which gives rise to 
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P 



(4.2) 
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With the help of these relations it is now easy to determine the remaining terms proportional to exp(ifcir) from all 
the possible combinations for q and p and the corresponding results for tpi calculated in appendix C. 

Using the scalar product ( p9| ) and the result for <po, the solvability condition at order 0(e 3 ) can be fo rmulated. 
It contains a term proportional to MiA^A^ which by eliminating Mi using the solvability condition (Cll) at order 
0(e 2 ) is transformed into terms proportional to A 2 | 2 ^4i and |A 3 | 2 A 1 . We then multiply the solvability condition at 
order 0(e 2 ) by e 2 and the one at order 0(e 3 ) by e 3 and add them together. Observing eM\ + e 2 M 2 — M — M c , 
returning to the original time by using e 2 d T = dt and introducing the scaled amplitudes A n = eA n we eventually end 
up with an amplitude equation of the form ([2.3£) with explicit expressions for the parameters j,gh,gt and g n . 



V. RESULTS 



The expressions for 7, gh, g t and g n are rather long and will not be displayed. Moreover, due to the large number of 
parameters in the two liquid system it is more appropriate to analyze some experimentally relevant parameter com- 
binations rather than to display cross sections along some direction of the parameter space. For the five experimental 
setups specified in appendix A the results of the non- linear analysis are summarized in the lower part of table 1. 

In order to finally address the planform selection problem we note that from the linear stability analysis of the roll, 
square and hexagon solutions of the amplitude equation (2.38|) it is well known p7],[l5f that: 



rolls are stable if gt > 1 , gt > 1 , gn > 1 , and e > jj 



"1 



(l-Sh) 2 ' 



squares are stable if 1 + g n < g^ + g t , \g n \ < 1, and e > ,^_J 



7"(l + gn) 



• hexagons are stable if 1 + 2gh > 0, e > en ■ 
1 + 2g h <g n + 2g t or e < e hts - —J-iSii+ML 



~ 



4(l+2 Sh )> 



9h-gt) 2 ' 
either g^ < 1 or e < ehtr 



n_ i2 1 and either 

l 1 9h ) 



(l+2g h -g n -2g t )i 



In addition to the special values of e defined in the last point above we have also included in table 1 the amplitude 
Ah of the pattern at onset. The he xago n pattern appears through a backward bifurcation which strictly speaking 
invalidates our perturbation ansatz ( 2.2E ). However, the interval of sub-critical hexagons as well as the amplitude of 
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the pattern at onset are for all investigated setups rather small such that the ansatz is still a good approximation for 
what really happens. 

Except for setup 3 when heated from below we always find g n < 1 excluding the possibility of stable rolls within 
the framework of our weakly non-linear analysis. For all setups we get 1 + 2gh > which implies that for hexagons 
the cubic term is able to stabilize the linear instability. Moreover for all setups gh > 1 and 1 + 2gh > g n + 2gt 
implying that the values of ehtr and ehts give the stability border for hexagons. Being the result of an expansion in 
the amplitude of the unstable modes the numerical values of ehtr and ehts are only reliable if they are not too large. 
If these values are hopelessly outside the validity of our perturbation approach they are not displayed in table 1. In 
all other cases we find ehtr > thts for setups with g n < 1 in accordance with the fact that rolls are then unstable to 
squares. The value of ehts is always positive which means that exactly at onset our analysis always predicts hexagons 
as the stable planform and excludes squares. However, in the cases where ehts is rather small (e.g. setup 4 when 
heated from below) hexagons get very quickly unstable to squares when passing the stability threshold. 

The sign of 7 is related to the detailed convection pattern of the hexagon planform. For 7 > the hexagons in 
the lower fluid are up-hexagons (liquid rises in the center) and the one in the upper layer are down-hexagons. For 
7 < the situation is reversed. We do not know of experimental results concerning this feature for the two liquid 
Marangoni problem. 





setup 1 


setup 2 


setup 3 


setup 4 


setup 5 


AT 


0.415 


4.032 


-3.945 


1.523 


-0.256 


0.859 


-18.957 


1.718 


rZ c 


2.495 


2.745 


0.714 


4.3416 


1.0328 


2.377 


0.861 


1.901 


M 


453 


1919 


-1878 


1978 


-333 


869 


-19188 


379 


R 


676 


654 


-640 


669 


-113 


733 


-16168 


45 


MW 


24.1 


614 


-601 


12107 


-2036 


149 


-3284 


777 


RW 


4.88 


592 


-579 


8145 


-1370 


49 


-1091 


143 


7 


0.406 


0.367 


-0.559 


-0.7478 


-0.5428 


0.423 


-0.507 


0.430 


9h 


1.225 


1.196 


1.411 


1.57 


1.36 


1.188 


1.417 


1.377 


9t 


1.442 


1.480 


1.501 


1.021 


1.529 


1.164 


1.273 


1.551 


9n 


0.030 


0.419 


0.075 


1.594 


-0.027 


-0.355 


-0.050 


0.628 


Qi 


-0.012 


-0.010 


-0.020 


-0.034 


-0.020 


-0.013 


-0.017 


-0.012 


A h 


0.118 


0.108 


0.1462 


0.180 


0.146 


0.125 


0.132 


0.114 


thtr 


- 


- 


6.30 


6.12 


7.50 


- 


5.04 


4.38 


this 


1.670 


- 


1.734 


7.922 


1.848 


0.180 


0.358 


- 



Table 1: Results for the critical temperature difference AT over both liquids (AT > for heating from below, AT < 
for heating from above), the critical wavenumber fc c , the Marangoni and Rayleigh numbers of both liquids at onset, 
the parameters of the amplitude equation ( 2.38| ), the sub-critical threshold e^ for the hexagonal pattern, its amplitude 
Ah at onset, and the values ehtr and ehts at which the hexagon pattern gets unstable towards the formation of rolls 
and squares respectively. If the numerical values of ehtr and ehts obtained are larger than 10 they are meaningless as 
result of a perturbation expansion in e and are therefore not displayed. 



For the parameters of setup 3 and a total depth of 4.5 mm we have again scanned the dependence of the results 
of the non-linear analysis on the thic kness of the bottom layer for the case of heating from below. FigJ6] shows the 
coefficients of the amplitude equation ( 2.38 ) as functions of hS 1 ' . The most apparent feature is the strong sensitivity of 
the coefficients on variations of the depth ratio. In experiments the depth must therefore be controlled very accurately 
in order to allow sensible comparison with the theory. The system under consideration shows a transition from up to 
down hexagons when varying the depth ratio as can be seen from the change of the sign of 7. 
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1.5 



0.5 



-0.5 





I 
1 






N 

N 


\ 
\ 
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- 


/' 
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- -j' \ 




\ 


\ 

i 


1 , 1 




i 



h 



(1) 



FIG. 6. The parameters gn (full), gt (dashed), g n (dotted) and 7 (dashed-dotted) of the amplitude equation ( 2.38 ) as 
functions of the thickness h, ' of the bottom layer for setup 3 with a total layer depth of 4.5 mm and heating from below. For 
1.5 mm< hr> < 2.5 mm the oscillatory instability precedes the static one. 

Finally in fig.fj] the dependence of Ehtr and £hts on /jW is displayed. For most values of h^ we have Ehtr > £hts 
and the hexagon pattern becomes unstable to the formation of squares. However for h^> = 1.5 also a secondary 
transition to rolls is possible. For some values of the depth ratio we find a very small £hts- Since at the same time 
also the absolute value of Sh is very small implying a small hysteretic window for the formation of hexagons it is quite 
conceivable that in these situations in the experiment the hexagon pattern cannot be observed at all and squares are 
seen directly at onset. 

1.5 



0.5 



-0.005 



-0.01 - 








h 



(1) 



FIG. 7. The values of Ehtr for the transition from hexagons to rolls (dotted) and Ehts for the transition from hexagons to 
squares (full) as functions of the thickness hs ' of the bottom layer for setup 3 with a total layer depth of 4.5 mm and heating 
from below. Also shown is the value Eh at which hexagons appear sub-critically. Note the different scales for positive and 
negative values at the vertical axis. For 1.5 mmm < tr 1 ' < 2.5 mm the oscillatory instability precedes the static one. 
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VI. DISCUSSION 

In the present paper a weakly non-linear analysis for Benard-Marangoni convection in systems of two superimposed 
liquids has been developed. A consistent treatment of the full hydrodynamics and heat conduction in both layers 
was performed. As crucial simplifying ingredient of our approach we have used the assumption of an undisturbed 
interface between the liquids. Comparison with the complete linear stability analysis including interface deflections 
has revealed that this approximation is extremely good for the pattern forming instability occurring at not too long 
wavelengths. We have considered the planform selection problem by determining the relative stabilities of roll, square 
and hexagon patterns. To this end the coefficients of the appropriate amplitude equation were calculated as functions 
of the hydrodynamic parameters by a perturbation theory in the amplitude of the unstable mode. As is well known 
|p0| , this expansion is not rigorous for the case of a sub-critical bifurcation leading to a finite amplitude immediately 
at onset. However, for the parameter combinations used we found that the hysteresis, as measured by e/,, is weak and 
the results obtained should therefore be rather accurate. 

Explicit numerical results were obtained for five different sets of experimentally relevant parameters of the fluids. 
Since the system is on the one hand characterized by nine dimensionless parameters whereas it is on the other hand 
very hard to find two really immiscible fluids to perform the experiments this seems to be the most sensible way to 
theoretically investigate the peculiarities of the system which may also be seen in experiments. For all parameter 
combinations investigated we predict hexagons at onset in agreement with recent experimental findings |22| . This 
shows that extrapolations from previous results on liquid-gas systems |f3] to the two liquid layer system which gave 
arguments in favour of squares directly at onset are potentially dangerous and the full hydrodynamics of both layers 
has to be taken into account. Moreover in most cases rolls were found to be unstable to squares for all values of the 
super-criticality parameter e. In particular for the parameter values of the experiments described in j2jj we do not 
find stable rolls in contrast to the secondary transition from squares to rolls reported for this case. 

The hexagonal pattern gets unstable to squares at a positive value ehts of the super-criticality parameter. For 
different experimental setups the values of thts differ substantially. Moreover even for the same combination of fluids 
it depends strongly on the depth ratio (cf. fig.0) ■ Nevertheless in most cases the values found are significantly smaller 
than those characteristic for Marangoni convection in single layer systems. In |32] the transition from hexagons to 
squares in an experiment with a single fluid layer were, e.g., reported to occur at e = 4.2 with the theoretical value 
resulting from a numerical integration of the Navier-Stokes equation being even higher. For the two-layer setups 4 
and 3 on the other hand studied in [M and p2| respectively ents is so small that it is well conceivable to miss the 
hexagonal pattern completely in the experiment and to observe squares as the first pattern after the instability in 
accordance with experimental findings. For setup 3 one also notes that together with ents also the absolute value of e/j 
characterizing the sub-critical stability region of the hexagon planform gets very small such that hexagons exist only 
in an extremely small window around criticality. Note also that our analysis is only concerned with perfect patterns 
hardly occurring in the experiments. It seems well possible that squares are generated by some inhomogencous 
nucleation process even before ehts is really reached. 

The remaining discrepancies between theoretical and experimental findings might be due to the perturbative char- 
acter of our derivation. In particular, there is the possibility of so-called asymmetric squares in pattern forming 
hydrodynamic systems [pT| which, bifurcating dis continuously from the quiescent state do not show up in a pertur- 
bative approachR. At the moment it is not clear whether these patterns can be expected already at the small values 
of the super-criticality parameter e used in the experiments. Since the flow pattern of asymmetric squares is rather 
different from the one of conventional squares it might be possible to clarify experimentally which form of squares has 
been observed. 
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would also like to thank F. Busse and W. Pesch for interesting discussions and Jean Bragard, Wayne Tokaruk and 
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Office of Life and and Microgravity Sciences Grant NAG3-1839. 
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APPENDIX A: PARAMETER VALUES 

In this appendix we have collected the values of the hydrodynamic parameters used for the numerical calculations 
of the present paper. All five sets correspond to experimentally relevant combinations. Setup 1-3 have been studied 
in |22] . Setup 4 was investigated in |21| whereas setup 5 is from the classical work 113] . 





setup 1 


setup 2 


setup 3 




lower fluid 


upper fluid 


lower fluid 


upper fluid 


lower fluid 


upper fluid 


substance 


HT135 


silicon oil 


HT 70 


silicon oil 


acetonitrile 


n-hexanc 


h(mm) 


2 


1 


.92 


2.14 


1.775 


2.725 


p(S) 


1730 


940 


1680 


920 


776 


655 


Kf) 


1 • 10~ 6 


2 • 10~ 6 


5-10- 7 


5 • 10~ 6 


4.76 -10- 7 


4.58 -10^ 7 


\ m.sK 1 


.070 


.134 


.070 


.117 


.188 


.120 


C P\TqK> 


962 


1498 


962 


1590 


2230 


2270 


<*(*) 


1.10 -to- 3 


1.10 -10" 3 


1.10- 10~ 3 


1.05 • 10- a 


1.41 -10- 3 


1.141- 10" a 


%{N/mK) 


-5-10~ 5 


-4.5 -lO" 5 


-l-io- 4 




setup 4 


setup 5 






lower fluid 


upper fluid 


lower fluid 


upper fluid 




substance 


FC75-FC104 


water 


water 


benzene 




h(mm) 


1.28 


2.78 


2.0 


1.0 




P&) 


1760 


998 


999 


885 




K^) 


8-10~ 7 


1 • 10~ 6 


1.14 -lO" 6 


7.77 -10" 7 




'"(msR-) 


.063 


.586 


.59 


.1615 




c p\~kW> 


1046 


4104 


4186 


1757 




a(-k) 


1.4 -lO -3 


2.06 -10~ 4 


1.50 -10" 4 


1.06- 10^ 




%{N/mK) 


-4.7 -lO" 5 


-5 • 10" 5 





Table 2: Parameter values for the five different experimental setups studied in this paper. In addition for all setups 
g = 9.81m/s 2 was used. Note that the value of da/dT is difficult to determine experimentally, the given values are 
therefore rough estimates or fitted from the linear analysis. 



APPENDIX B: OPERATOR EXPANSION AND ADJOINT PROBLEM 



The decomposition (2.32) of the linear operator is not completely straightforward for the Marangoni problem 
because the bifurcation parameter M not only occurs in the linear operator but also in the corresponding boundary 
conditions. A transparent way to deal with the situation is to include the boundary condition involving M into the 
operator L |7J], which is then written in the form 



L = 



A 2 
1 





cMA± 
A 










vA 2 
l 



-vd%=o 







acMAj 

X A 




acting now on the correspondingly augmented state vector 



ip = 



( W \ 
W 

e 



\ 



o 





-maJ 



(Bl) 



(B2) 
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The operator is completed by the boundary conditions 

w = d z w = 6 = at z = — 1 
w = W = , d z w = d z W ,e = Q,d z = nd z e . 
W = d z W = 6 = at z = a 



at 



(B3) 
(B4) 

(B5) 



which differ from ([2.23])-(p.25]) just by the omission of the boundary condition involving M. We now easily find 



L 



and 



/ A 2 


cM c A± 










o \ 


1 


A 
















^A 2 


acM c A± 








l 




X A 


V£l»=o 





-Vd% 


=0 


-M c Aj 


/<> 


cMiA_l 








° \ 





















Li = 








acMiA_L 





) 





















Vo 











-M X AJ 


/° 


cM 2 A_l 








° \ 





















L 2 = 








acM 2 A± 





) 























Vo 











-m 2 aJ 





(B6) 



(B7) 



(B8) 



where all three operators are completed by the boundary conditions (B3)- 

The adjoint operator is defined by ((p\Lip) = (L + ip\ip). Introducing the scalar product 



1 



,L/2 ,L/2 

(<P\<P) = r lim T2 I dx I d y 

L^oo L s J_ L/2 J-L/2 



dz{w*w + 6*6) 



dz{W*W + 9*9) + d z w*\ z = o 



z=0P 2=0 



(B9) 



we find after some partial integration that L + is given by 



L+ = 



A 2 


1 








\ 


cMAj_ 


A 

















^A 2 


J_ 











acMA ± 


X A 





o 


—9 z \ z= o 





xd z \ z =o 


-MAJ 



(BIO) 



acting on the augmented vector 



/ 



w 



<P 



\ 



and completed by the boundary conditions 

w = d z w = 

W = 0,d z w 



at z 
1 



W 

e 

\d z w\ z=0 J 



1 



(Bll) 



w 



X, 



d z W , dlw = vdlW , 6=^Q. 

p K 

W = d z W = 6 = at z = a . 



at z = 



(B12) 
(B13) 
(B14) 



It is, of course, possible to transform back the last line of L + into a boundary condition and this is indeed advantageous 
to determine (p~o explicitly, however for the use in the solvability conditions the above augmented form is the most 
appropriate one. 
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APPENDIX C: THE 0(e 2 )-PROBLEM 



In this appendix we solve eq.( |2.34[ ) for the case of a static instability. From the term JV(ipo, fo) and the structure 
(2.37) of ipo it is clear that the right hand side of this equation will contain several terms with different exponential 
factors of the form exp(i(±fe n ± k m )r). Because of the linearity of the equation we may solve it separately for all 
these term in the inhomogeneity. 

Let us start with the so-called non-resonant t erm s in which the angle <fi between ±k n and ±fc m is differ ent from 
2tt/3. It is clear then from the z-y-integrals in (B9) that for these terms {(po\J\f(ipo, (po)) = 0. In view of ( ] B7|) the 
solvability cond ition boils down to M\ = and hence removes the Li^o-term from the inhomogeneity of ( |2.34|) . Using 
the form ([2.37|) of tpo we therefore find as equations for ip\: 



A 2 w 1 +cM c A ± 6 1 =A n A m e l( - ±k " ±k ^ r — [(1 + cos(0))(u#'wo + (1 - 2cos($K«#) - 2£; 2 sin 2 ^) u, ^)] 
wx + A0! = A n A m e l{ - ±k - ±k ^ r 2 (w 9' - cos(0) w' 6 ) 
vA 2 W x + acMAxQi = A„A m e l(±fc " ±fcm)r J- [(1 + cos(<£))(W^"Wo + (1 - 2 cos((f>))W^') - 2fc c 2 sin 2 (</>) W W&) 



-Wi + xA6! = A n A m e l( - ±k - ±k ^ r 2 (W & - cos(0) Wfoo) 



where the prime denotes d iffere n tiatio n with respect to z. Since M\ — the boundary conditions completing this set 
of equations are given by ( 2.23 )-( 2.25 ) with M — M c . 

The solution of these equations is of the form </?i = A n A m (fi(z) exp(i(±fc„ ± k m )r). We first determine a solution 
of the inhomogeneous equations using the ansatze 



wT\z) 



6 

E 



wuje 



(Xi+Xjjz 



12 

W p h (z) = £ w Uj e 



(A t + A,)z 



6™\ Z ) = J2 9 Uj e 



{\i+\j)z 



12 



e\ nh (z) = J2 hiJ e 

i,j=7 



(A t + A,)z 



(CI) 
(C2) 



which give rise to algebraic equations for the coefficients Wuj, Ouj, Wuj and @uj in terms of woi and Aj. This solution 
does not yet satisfy the boundary conditions. We therefore add a proper solution of the homogeneous equation which 
is written in the form 



w^ om (z) 
W? om (z) 



\ * -..horn „\iZ 






W 



It 



1 = 1 
12 

i=7 



qhorn 



(z) = - 



6 



-A 2 -2fc2(l + cos(^)) 



ei lom (z) 



i ^ 



..horn 



,A,z 



«X^Af-2fc c 2 (l + cos(0)) 



(C3) 

(C4) 



with Xi satisfying 



[A 2 -2fc 2 (l + cos(0))] 3 = 



-2cM c fc 2 (l + cos(</>)) for i = 1, ... ,6 

-2-f-cM c fc c 2 (l + cos(0)) for $ = 7, . . . , 12 



(C5) 



Note that A^ ^ Xi . Therefore the determinant of the matrix in the inhomogeneous set of linear equations for u^° m 
is different from zero and the solution is unique. Note also that for <j> = n the procedure can be simplified since 
wi(z)=Wi(z) = 0. 

As for the resonant terms arising from the interaction of modes with an angle <fi = 27r/3 between their respective 
±fc-vectors let us focus on the one proportional to exp(ifcir). It is has one contrib ution proportional to A\ stem ming 
from —Litpo and another one proportional to A^A^ originating from AA(^o, <Po) in (2.34). Using L\ as defined by (B7) 
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the resulting equations are of the form 



AV + cM c Aj.(9i = e lkir ^^-(w "w + 2w'w" - 3k 2 c w Q w') + A x cM x k 2 6 

Pr 

w x + A0i = e ihir A* 2 A% (2 w 9' + w' 8 ) 



vt\ 2 W x + acM A±Q X 



A* A* 

Jkxr A 2 A Z 

Pr 



{Wg'Wo + 2 W^Wg - 3k 2 W W^) + A x acM x k 2 9 



1. 



W x + X A6i = e lkir A*A* 3 (2 W Q Q' + W^@ ) 



(C6) 
(C7) 

(C8) 

(C9) 



The boundary con diti ons are again given by ( p.23| )-( [2.25| ) except for the one containing the Marangoni number, which 
is modified to (cf. B7) 



d 2 w x - r]d 2 W x - M C A ± 6 X = -A x e lkir M x k% at 







(CIO) 



Due t o th e resonant factor e lkir the terms arising from Af((po,(po) are not automatically perpendicular to (po and 
using (B9) the solvability condition acquires the non-trivial form 



= A*A* 3 



+A X cM x k 2 



w 
<!: ( -p^K'^o + 2w' w'o ~ 3k 2 c w w' )+6* (2w e' + w' e ) 

w* 

dz [ -J^( W o"Wo + 2 W^Wg - 3k\ WWo) + 65(2 W Q' + W^Q ) 



dzwQ0 o + a / (IzWqQq d z Wo\ z=Q 6 \ z=0 



(Cll) 



We use this equation to replace the terms involving M x in eqs.(|C6|)-(|C9|) and in the boundary condition ( CIO ). The 
solutions to these equations can then be written in the form ^4^3 ViO*) e lklV . Again we first determine a particular 
solution of the inhomogeneous equations by using the ansatze: 



< nh (z) = E w ^ e(A,+ " ,)2 

12 

Wl nh {z) = J2 w XlJ e^ +x ^ z 

i,3=7 



E 

(=i 

12 

E 

i=7 



,A» 



w X i ze 



,A, 



w X i z e 



n nh (z) = E ^ (Ai+A ^ + E( 01 * z + *») e> 

i,j=l i=l 

12 12 

hi 3 e { 

i,3 = 7 



e\ nh (z) = J2 d XlJ e^ +x ^ z + J2(0u z + 0u) e Al 



i=l 



To s a tisfy the boundary conditions we add a solution of the homogeneous equations which must be of the form (cf. 

©,©) 



lt hoin i 



horn p\iZ 



;) = E-i 

i=l 
12 



qhom i 



e^ om (z) 



2-^ \2 _ 



.horn 



,\iZ 



4=1 A * 



k 2 
horn 



i=7 



1 12 

1 y- ^i. 

KV 2-^1 \2 _ 1,2 



(C12) 
(C13) 



The boundary conditions give rise to an mhomogeneous system of linear equations for the coefficients w^° m with 
the same singular matrix A which appeared in the linear stability analysis. Due to the solvability condition ( Cll ) 
however, the inhomogeneity of this set of linear equations is perpendicular to the zero eigenvector of the adjoint 
problem and therefore the system admits solutions. Their numerical determination is most conveniently done by 
using the singular value decomposition of the matrix A p9| . This method yields an approximate solution even if 
the solvability condition is not fulfilled exactly, which will always be the case due to numerical errors. Moreover, 
the so-called residual quantifying the deviation from the exactly solvable case gives another check of the numerical 
accuracy of the whole procedure. 

Finally, the solution for w x ° m obtained in this way is not unique since one can always add a solution of the 
homogeneous equations. We will enforce the additional constraint 



0= (<A)|v?i) 



lim —r 

L^oo L 2 



L/2 pL/2 

dx I dy 

-L/2 J-L/2 



dz( w * w x + e* e x ) + / dz{W£w x + e* e x ) 

! JO 



(C14) 
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to remove this ambiguity. The rationale behind this re quirem ent is as fo llows. Assume that we knew the exact 



quircm 



solution ip of the full non- linear problem. According to ( |2.2S| ) and ( 2.37 ) we want A n to be the a mplit ude of the 



contribution to <p proportional to exp(ik n r), i.e. (exp(ik n r)(po(z)\ip) = eA n . Using the e xpa nsion ([2.29) for ip this 



results in ((po\tpi) = for all Z > 1. Note the use of different scalar products in (C14) and (B9). 

This completes the solution of the 0(e 2 ) equations. The results are specified by the various matrices 
wuj,w u , <°"\ BujJu, hi and 6^ m - 
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